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Emergence of coexisting percolating clusters in networks 
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It is commonly assumed in percolation theories that at most one percolating cluster can exist in a network. We 
introduce sausage-like networks (SLNs), an ensemble of synthetic modular networks in which more than one 
percolating cluster can appear. We show that coexisting percolating clusters (CPCs) emerge in such networks 
due to limited mixing, i.e., a small number of interlinks between pairs of modules. We develop an approach 
called modular message passing (MMP) to describe and verify these observations. We demonstrate that the 
appearance of CPCs is an important source of inaccuracy in the previously introduced percolation theories, such 
as the message passing (MP) approach. Moreover, we show that the MMP theory improves significantly over 
the predictions of MP for percolation on synthetic networks with limited mixing and also on several real-world 
networks. These findings have important implications for understanding the robustness of networks and in 
quantifying epidemic outbreaks in the susceptible-infected-recovered (SIR) model of disease spread. 
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Percolation theories are among the most studied in net¬ 
work science [1], as well as in several other areas [2, 3], pro¬ 
viding insights for a broad range of applications such as ro¬ 
bustness of a network to random failures or attacks [4], epi¬ 
demics in contact processes [5], vaccination strategies [1], 
neuronal avalanches [6], and stability of gene regulatory net¬ 
works [7]. In the simplest case of bond (or site) percolation, 
a fraction p of the links (nodes) are randomly chosen to be 
occupied and the rest of the links (nodes) are removed from 
the network [1]. The quantity of interest is S: the expected 
fractional size of the giant component (GC) of the network, 
which in the limit of infinitely large networks is referred to 
as the percolating cluster (PC) of the network [2, 8], The 
size of the GC scales linearly with the network size while 
the fractional sizes of other clusters vanish in the limit of 
infinitely large networks. 

Theoretical approaches and extensive numerical simula¬ 
tions play pivotal roles in understanding and describing the 
behavior of percolation processes on networks. The pi^ the¬ 
ory for bond percolation [1, 9], for example, can accurately 
describe the results of numerical simulations on configura¬ 
tion model [10] networks using only the network degree dis¬ 
tribution. On networks with degree-degree correlations, the 
accurate results are obtained using the so called P{k, k') the¬ 
ory [11] which employs the joint degree distribution. The 
theory [12] can provide a more accurate description of 
dynamics on modular networks, as it considers the joint de¬ 
gree distributions within and between modules. The mes¬ 
sage passing (MP) approach [13] provides more accurate 
results than the aforementioned theories as it uses the full 
information on the adjacency of individual nodes, and re¬ 
duces to the above degree-based approximations in special 
cases [14]. 

As these theories assume the network is locally treelike, 
they are prone to errors in clustered networks which have 
an appreciable density of short loops [9, 14]. However, on 
some real-world clustered networks these theories still per¬ 
form well, and in some other cases the inaccurate predictions 
of these theories are shown to be only partly caused by the 
presence of short loops [14]. This indicates the presence of 
an unexplained source of error and possibly a phenomenon 


not captured by the theories. 

In this paper, we show the appearance of coexisting perco¬ 
lating clusters (CPCs) in certain networks, and demonstrate 
that this phenomenon causes significant errors in the afore¬ 
mentioned theories. We show that CPCs appear in modular 
networks with limited mixing, i.e., networks with a suffi¬ 
ciently small (limited) number of interlinks between mod¬ 
ules. We verify these observations by developing the modu¬ 
lar message passing (MMP) theory which takes into account 
the presence of independent CPCs. We show that the MMP 
theory provides accurate predictions on treelike modular net¬ 
works with limited mixing and also improves over the pre¬ 
dictions of MP on several real-world clustered networks. 

We begin by introducing sausage-like networks (SLNs), 
a simple ensemble of random networks that demonstrates 
the appearance of CPCs. To create an SLN we first pick 
a graph with size N^, which can be any connected undi¬ 
rected unweighted graph. Then we make M identical copies 
of that graph (Fig. la) which will become modules in the 
SLN. We assign to each of these modules a unique label 
m 6 (1,2,... ,M] and connect each pair of modules with 
consecutive labels m and m -H 1 by / links. To do so, exactly 
7/2 links are selected randomly from module m. For each se¬ 
lected link im — jm we consider its copy — jm+i in mod¬ 
ule m+\, and rewire these two links to create two new links 
i,„ — i,„+i and —y„, instead. The resulting SLN is com¬ 
prised of a chain of modules (Fig. lb), each pair of consecu¬ 
tive modules connected with exactly I interlinks. Moreover, 
an SLN has a degree distribution and degree-degree corre¬ 
lations between and beyond the nearest neighbors identical 
to those of the original graph [14]. Similarly, one constructs 
SLNs from non-identical modules by rewiring links that are 
not copies but randomly selected from each module. 

In Figs. Ic-le, we illustrate the bond percolation results 
on SLNs constructed from a 3-regular graph. The results 
from the MP theory are obtained using the Eqs. (1) and (2) 
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FIG. 1. (Color online) (a) A schematic of M = 3 identical copies of 
a graph, (b) that are connected in a chain to construct an SLN with 
M modules, each pair of modules connected through 1 interlinks, 
(c) Bond percolation results on SLNs constructed from a 3-regular 
graph of size Nm = 1000 for different /. (d) The mean absolute 
error E versus I and (e) versus N„ for the MP theory on such 3- 
regular SLNs. In (c) M = 50, in (d) = 5000, and in (e) M = 10. 


of Ref. [14] (originally due to [13]): 


Uij ^ I - p + pY\AjkUjk, (1) 

k^i 



where p is the occupation probability, is the probability 
that node i is not connected to the network PC via its link 
to j, A is the network adjacency matrix, N is the size of the 
network, and i, = 1 “ H / the probability that node i 

is in the PC. 

Numerical simulations show that the behavior of S in the 
SLNs depends substantially on the number of interlinks I 
and the number of modules M. Surprisingly, the result of 
the MP theory for the SLNs is independent of I and M and 
coincides with its prediction for a single 3-regular graph 
(Fig. Ic); it is worth mentioning that in such SLNs the result 
of MP theory is the same as the results of any of its degree- 
based reductions mentioned above. The numerical results, 
on the other hand, deviate from the theoretical prediction as 
the number of interlinks I is decreased (Fig. Ic). We quan¬ 
tify the difference between the theoretical and numerical re¬ 
sults by calculating the mean absolute error between the two: 
E = \Stheo{Pj) - ~Snum{Pj)\, where the sum is over 

R - 100 equally spaced occupation probabilities pj - HR. 

Figure Id shows that the error increases dramatically for 
lower number of interlinks /; on the other hand the error 
increases only slightly for smaller module sizes N,n (Fig. le). 
Figure le highlights also that the error increases for lower I 
rather than for lower ratio of I to all edges. 

To understand these observations, it is necessary to in¬ 
spect also the numerical results for S , the fractional size of 
the network PC for single realizations of the bond perco¬ 
lation process. Figure 2a shows the results for a 3-regular 
SLN with M - 5 modules; the MP theory overestimates S, 


which is the expected (average) value of S. For a sufficiently 
large network, we normally expect that S fluctuates slightly 
around S. Surprisingly, we observe in Fig. 2a that, for a fixed 
value of p (in a certain interval), S can take one of several 
possible values which can be significantly different from S. 
The different possible values for S can be explained as fol¬ 
lows. Let us denote by 5 m the fraction of nodes in module 
m that are in the network PC predicted by the theory. Hence, 
for the MP theory Yiieih where m denotes the set 

of nodes located in module m, and Eq. (2) becomes: 

M 

(3) 

m=l 

In Fig. 2a, we can observe that S predicted by MP matches 
the largest possible value of S. The lowest possible value of 
5, on the other hand, coincides with Sm of only one module. 
The next 3 larger possible values of S coincide with the sum 
of 5m for 2, 3, and 4 modules respectively (Fig. 2a). 

Consider another example of an SLN consisting of one 
3-regular module and one 4-regular module each having the 
same number of nodes Nm, and let us denote by 5 1 and S 2 
the 5 m of the 4-regular and 3-regular graph respectively. For 
this SLN, 5 from numerics and from the MP theory both 
match 5 1 up to p-ireg - 0.5, the percolation threshold of a 
3-regular graph. Above this value, the MP prediction devi¬ 
ates from the numerical result (Fig. 2b). This deviation can 
be better understood by looking at single realizations of the 
Newman-Ziff algorithm [ 8 ], where starting with no occupied 
links, we occupy links one by one in random order. As the 
network is large, p is approximately equal to the fraction of 
occupied links. In Fig. 2b, we can observe that up to p^reg^ 
single realization values of 5 match the value of 5. However, 
above p^reg, while the predicted 5 is 5 1 H- 52 , 5 will remain 
equal to 5 1 until some larger value of p and then suddenly 
jump to 5 1 H- 52 . 

This implies that although 5 2 is finite for p > p^reg, it does 
not represent the probability a node in module 2 belongs to 
the network PC, and accordingly 5 1 H- 5 2 is not the expected 
size 5 of the network PC. In fact, the nodes in module 2 are 
part of a PC with size S 2 and those in module 1 are part 
of another PC with size 5 1 . In a single realization these 
two PCs may be unconnected. Then, when more links are 
occupied one at a time, at a value of p with p > Pireg the 
two PCs become suddenly connected, and the size of the 
network PC changes abruptly from 5 1 to 5 1 - 1 - 52 . We refer 
to such PCs, which exist independently of each other in the 
network but may only be connected with a finite probability, 
as coexisting percolating clusters (CPCs). 

Next, we develop the modular message passing (MMP) 
theory to describe and verify the phenomenon of CPCs ap¬ 
pearing in networks. The two main assumptions are (i) mod¬ 
ules can percolate independently (hence, the appearance of 
independent CPCs), and (ii) PCs of neighboring modules 
are connected with probability 77 < 1. This is a new con¬ 
cept different from the common assumption (see for exam¬ 
ple [1, 12, 15]) that there exists only one monolithic PC in a 
network; here we show that in networks with limited mixing, 
the network PC is polylithic, i.e., constituted by CPCs that 
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FIG. 2. (Color online) (a) Percolation on an SLN with M - 5 and 
= 20000. For each p, S is shown separately for 20 single re¬ 
alizations of percolation (triangles). S for single realizations does 
not match S (circles) averaged over 500 realizations, nevertheless it 
coinsides with one of the 5 possible values denoted by the dashed 
line (MP theory) or the solid lines, (b) Percolation on an SLN 
with a 3-regular and a 4-regular module, 1-4, and N„, = 20000. 
Blue lines represent the results from 20 single realizations of the 
Newman-Ziff algorithm. The vertical line on the left (right) de¬ 
notes the percolation threshold of a 4-regular (3-regular) graph, (c) 
The MMP prediction for S on the SLN of panel (b) matches the 
numerical results on this SLN (circles) as well as on similar SLNs 
with identical I and any sufficiently large N„,. (d) Predictions for 
the SLN of panel (a). The 5th and 95th percentiles of S are in¬ 
dicated by the green shade (numerical results) and by the dashed 
lines (MMP theory). 


are connected together. In such networks, the CPCs emerge 
independently inside the internally well connected modules, 
but they still may not be connected to each other due to the 
small number of interlinks between the modules. For p < I, 
interlinks may be unoccupied; the CPCs of two neighboring 
modules are connected if and only if they share at least one 
occupied interlink that connects nodes from the CPCs of the 
two modules. 

For SLNs we assume that each boundary node (a node 
with links to other modules) has exactly one link to a neigh¬ 
boring module. Then for an SLN with two modules, our 
MMP theory is comprised of two simple equations. First 
we calculate the probability that the CPCs of the two 
modules m and n are connected: 


rjinn ^ I - (I - PVmVn)'■ (4) 

Here p is the occupation probability and is the probability 
that a boundary node in m is part of the CPC of m. Then 1 - 
pVmVn is the probability that the two CPCs are not connected 
via an interlink; raising this term to the power of I gives 
the probability that they are not connected via any of the I 
interlinks. Equation (4) is independent of Nm, hence in the 
thermodynamic limit {N oo), if I is fixed to a sufficiently 
small number, (independent) CPCs emerge, connected with 
a probability rjmnip) < 1 ■ Whereas for large /, since rjmn 1 
the two CPCs are connected with high probability, leading to 
a monolithic PC. For the SLNs described above, and v„ 
can be simply calculated using the p^ theory (see Sec. S.l. 


of the Supplemental Material (S.M.) [16]). For an SLN with 
two modules, the size of network PC is then: 

= ?7i2 (‘^1 +‘^ 2 ) + (1 -'7i2)‘5i> (5) 


where S \ and S 2 are the fractional sizes of respectively the 
larger and the smaller CPC of the SLN. Figure 2c shows 
that the prediction of the MMP theory (Eqs. (4)-(5)) matches 
perfectly the numerical result for the SLN of Eig. 2b. 

Eor SLNs with more than two modules Eq. (4) can still 
be used to calculate rj^n for each pair of modules m and n. 
In the case when boundary nodes have more than one inter¬ 
link or when the connection pattern of the modules can not 
be well approximated using the theory, Eq. (4) should be 
extended to include more information on the network struc¬ 
ture. We can use the full information on the adjacency of 
individual nodes to write a general formula for the connec¬ 
tion probabilities r]mn between CPCs in a treelike network: 





\ 




\ 

1 I I U[j( 

, kej{„(i) J 


(6) 


where i is a boundary node of module m and denotes 

the set of neighbors of i in modules m. Here, Uij and 
are the probabilities that i is not connected to the CPC of, 
respectively, modules m and n via its links to nodes inside 
each of those modules. Hence, within the hrst (second) set 
of parentheses in Eq. ( 6 ) is the probability that i is in the 
CPC of m (the CPC of n), and in the square brackets we 
have the probability that the two CPCs are connected via the 
interlinks of i. Therefore, is the probability that the two 
CPCs are connected via any of their I interlinks. 

To calculate S for networks with more than two modules 
(and consequently more than two CPCs), Eq. (5) should be 
extended as well. If a networks contains several CPCs then 
in a single realization of percolation, different CPCs (with 
sizes S,n) are connected together with some probability, cre¬ 
ating larger polylithic PCs. The polylithic PC I has a size 
^(poi) _ 5 ^^ where I denotes the set of CPCs that con¬ 

stitute 1. Then in a single realization at a hxed value of p, the 
size S of the network PC is max/ ^5 i.e., the size of the 

largest polylithic PC. Then the expected size of the network 
PC is 


S^YuP^s)S, (7) 

s 

where P^s) is the probability that in a single realization the 
size of the largest polylithic PC is 5. To calculate P(s), for 
each p we hrst calculate S^ip) values using the MP theory 
or an appropriate degree-based reduction of MP. Then we 
assume a meta-network which is comprised of meta-nodes; 
each meta-node m represents a module of the original net¬ 
work and has a weight S„,{p)- In each realization of the 
meta-network for a hxed value of p, each pair of meta-nodes 
m and n are connected with probability rjmnip). Eor suffi¬ 
ciently large number of meta-network realizations, we cal¬ 
culate P(j) using the Newman-Ziff algorithm [ 8 ] with the 
following modihcations: (i) each link is added with prob¬ 
ability rjmnip)^ and (ii) the size of a cluster I, comprised of 
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FIG. 3. (Color online) (a) Bond percolation results on an LFR net¬ 
work with N nodes, degree distribution P(k) oc community 
size distribution P(N„) oc N^', average degree k^^, the mixing pa¬ 
rameter ji, and minimum and maximum community sizes of re¬ 
spectively 50 and 1000. (b) Results for the western US power grid 
network. For the LFR network, MMP predicts accurately S and 
also percentiles for the distribution of S ; for the power grid its 
predictions match the results for the rewired version of this 
network. 

connected meta-nodes, is the sum of the weights of meta¬ 
nodes it includes, i.e., These calculations are per¬ 

formed very quickly, as the number of the modules is usually 
much smaller than the number of nodes. 

Figure 2d illustrates that, using Eqs. (6)-(7), the MMP the¬ 
ory performs very well for the SLN of Fig. 2a. As shown 
in Figs. 2a and 2b in the presence of CPCs the values of 
S for single realizations can deviate considerably from the 
expected value S. In the MMP theory, the variability of S 
is determined by defined above. The percentiles cal¬ 
culated from P(s) match very well the numerical values in 
SLNs (see Fig. 2d for example), which confirms that the 
high variability of S originates from the presence of CPCs 
that suddenly merge together. Similar results (not shown) 
are obtained for SLNs with different values of /, M, Nm, and 
also for SLNs constructed from modules with a heteroge¬ 
nous structure, i.e., with a power-law degree distribution. 

We also provide results for LFR benchmark networks 
[17], as an example of modular networks with heteroge¬ 
neous structure. In LFR networks, all pairs of modules can 
be connected without any restriction (as apposed to SLNs), 
and the node degrees and community (module) sizes have 
a power-law distribution. On LFR networks with low mix¬ 
ing between modules, the numerical results do not match 
the MP predictions, whereas MMP method provides accu¬ 
rate predictions (e.g., see Fig. 3a). As mentioned before, 
MMP assumes that modules have (independent) CPCs. This 
assumption holds only for sufficiently large modules. Never¬ 
theless, if the size of the modules is as small as ^ 50, 
MMP still performs well (see S.M., Sec. S.2. [16]). 

We can show that on several real-world networks the bond 
percolation results are affected by the emergence of CPCs. 
To do so, we first identify the best representations of net¬ 
work modular structure using a multiresolution community 
detection method (e.g., see Refs. [18-21]). Then, we choose 
the representation that maximises the modularity Q [22] and 
minimizes the number of modules with sizes smaller 
than = 50, i.e., the representation with maximum 

Q. Figure 3b shows that MMP improves significantly 
over the MP prediction of S on the western United States 
power grid network [9, 23], and also provides a prediction 


for the variability of S according to the percentiles of P(s). 
The results of MMP on the power grid network, match the 
numerics for the rewired [12] version of this network 
in which the links are rewired inside each module; 
rewiring preserves the modular structure and degree-degree 
correlations but effectively destroys the short loops. This 
shows that MMP provides a highly accurate prediction in 
the absence of short loops and when the modular structure is 
identified accurately. Similar results are shown in S.M. for 
several other examples of real-world networks [16]. 

In summary, we demonstrated that CPCs can emerge in 
ensembles of random networks and in real-world networks, 
when the network modules are connected via a small number 
of interlinks. Moreover, we showed that CPCs are an impor¬ 
tant source of error in the theories for bond percolation and 
proposed the MMP theory that accurately captures the im¬ 
pact of CPCs on percolation results. An important implica¬ 
tion of the appearance of CPCs is the uncertainty they cause 
in determining the network robustness: when the CPCs 
emerge, the size S of the network PC can be highly variable 
(as shown by the percentiles). This implies the prominent 
role of interlinks in network robustness, even in the absence 
of module-based targeted attacks [24], and subject to only 
random failures. Another implication is that the eventual 
size of an epidemic spread in the SIR model [25, 26] may 
not be best represented by S, which is the expected size of 
the largest polylithic PC, since even CPCs not in the largest 
polylithic PC represent (independent) outbreaks of compara¬ 
ble sizes located in different modules. Hence, the total size 
of an epidemic outbreak may better be represented by the 
sum of the sizes of all CPCs. 

We would like to thank Marian Boguna, Mel Devine, 
and David O’Sullivan for valuable comments. We acknowl¬ 
edge funding from Science Foundation Ireland (under pro¬ 
grammes ll/PI/1026 and 12/IA/1683) and from the Eu¬ 
ropean Commission EET-Proactive project PLEXMATH 
(EP7-ICT-2011-8; Grant No. 317614). 


[1] M. E. J. Newman, Networks: An Introduction (Oxford Uni¬ 
versity Press, Oxford, 2010). 

[2] K. Christensen and N. R. Moloney, Complexity and criticality 
(Imperial College Press, 2005). 

[3] D. Stauffer and A. Aharony, Introduction to percolation the¬ 
ory (CRC press, 1994). 

[4] D. S. Callaway, M. E. J. Newman, S. H. Strogatz, and D. 1. 
Watts, Phys. Rev. Lett. 85 , 5468 (2000). 

[5] M. A. Munoz, R. Juhasz, C. Castellano, and G. Odor, Phys. 
Rev. Lett. 105 , 128701 (2010). 

[6] E. J. Friedman and A. S. Landsberg, Chaos 23 , 013135 
(2013). 

[7] S. Squires, E. Ott, and M. Girvan, Phys. Rev. Lett. 109, 
085701 (2012). 

[8] M. E. J. Newman and R. M. Ziff, Phys. Rev. E 64, 016706 
( 2001 ). 

[9] S. Melnik, A. Hackett, M. A. Porter, P. J. Mucha, and J. P. 
Gleeson, Phys. Rev. E 83 , 036112 (2011). 

[10] M. Molloy and B. Reed, Random Structures and Algorithms 
6, 161 (1995). 



O Western US power gird network 
-MP 
-MMP 

— 95th and 5th percentiles 

□ P'l‘,, rewired network 
kk ' 

5th-95th percentile ✓ 













5 


[11] A. Vazquez and Y. Moreno, Phys. Rev. E 67, 015101(R) 
(2003). 

[12] S. Melnik, M. A. Porter, P. J. Mucha, and J. P. Gleeson, Chaos 
24 , 023106 (2014). 

[13] B. Karrer, M. E. J. Newman, and L. Zdeborova, Phys. Rev. 
Lett. 113 , 208702 (2014). 

[14] A. Eaqeeh, S. Melnik, and J. P. Gleeson, Phys. Rev. E 91 , 
052807 (2015). 

[15] P. Colomer-de-Simon and M. Boguna, Phys. Rev. X 4 , 041020 
(2014). 

[16] “See Supplemental Material at https://goo.gl/w5uwwo for 
some details on the calculations and additional examples for 
the MMP results.” . 

[17] A. Lancichinetti, S. Eortunato, and E. Radicchi, Phys. Rev. E 
78,046110(2008). 

[18] P. Ronhovde and Z. Nussinov, Phys. Rev. E 80 , 016109 


(2009). 

[19] V. D. Blondel, J.-L. Guillaume, R. Lambiotte, and E. Lefeb- 
vre, J. Stat. Mech. 2008 , P10008 (2008). 

[20] A. Eaqeeh and K. A. Samani, Physica A 391, 2463 (2012). 

[21] M. Rosvall and C. T. Bergstrom, Proc. Nat. Acad. Sci. U.S.A. 
105 , 1118 (2008). 

[22] M. E. J. Newman and M. Girvan, Phys. Rev. E 69 , 026113 
(2004). 

[23] D. J. Watts and S. H. Strogatz, Nature (London) 393 , 440 
(1998). 

[24] B. R. da Cunha, J. C. Gonzalez-Avella, and S. Gonsalves, 
arXiv: 1502.00353 (2015). 

[25] M. E. J. Newman, Phys. Rev. E 66, 016128 (2002). 

[26] B. Karrer and M. E. J. Newman, Phys. Rev. E 82, 016101 

( 2010 ). 


